﻿using System;
using System.Data;
using System.Data.SqlClient;
using System.Data.SqlTypes;
using Microsoft.SqlServer.Server;
using SmartMathLibrary.SQLServer.Commands;

public partial class UserDefinedFunctions
{
    /// <summary>
    /// The standard integration value for the special functions.
    /// </summary>
    private const int StandardIntegrationValue = 1000000;

    /// <summary>
    /// The Runge function.
    /// </summary>
    /// <param name="x">The specified x position for the function.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double RungeFunction(double x)
    {
        return 1/(Math.Pow(5*x, 2));
    }

    /// <summary>
    /// The Runge function. Alias command for <see cref="RungeFunction"/>.
    /// </summary>
    /// <param name="x">The specified x position for the function.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double Runge(double x)
    {
        return UserDefinedFunctions.RungeFunction(x);
    }

    /// <summary>
    /// The Heaviside function.
    /// </summary>
    /// <param name="x">The specified x position for the function.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static int HeavisideFunction(double x)
    {
        if (x <= 0)
        {
            return 0;
        }

        return 1;
    }

    /// <summary>
    /// The Heaviside function. Alias command for <see cref="HeavisideFunction"/>.
    /// </summary>
    /// <param name="x">The specified x position for the function.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static int Heaviside(double x)
    {
        return HeavisideFunction(x);
    }

    /// <summary>
    /// The Gudermann function.
    /// </summary>
    /// <param name="x">The specified x position for the function.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double GudermannFunction(double x)
    {
        return 2*Math.Atan(Math.Tanh(x/2));
    }

    /// <summary>
    /// The Gudermann function. Alias command for <see cref="GudermannFunction"/>.
    /// </summary>
    /// <param name="x">The specified x position for the function.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double Gudermann(double x)
    {
        return GudermannFunction(x);
    }

    /// <summary>
    /// The Lorentzcurve function.
    /// </summary>
    /// <param name="x">The x value for the function.</param>
    /// <param name="t">The centrum t.</param>
    /// <param name="s">The breadthsparameter s.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double LorentzCurveFunction(double x, double t, double s)
    {
        return (1/Math.PI)*(s/(Math.Pow(s, 2) + (Math.Pow(x - t, 2))));
    }

    /// <summary>
    /// The Lorentzcurve function. Alias command for <see cref="LorentzCurveFunction"/>.
    /// </summary>
    /// <param name="x">The x value for the function.</param>
    /// <param name="t">The centrum t.</param>
    /// <param name="s">The breadthsparameter s.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double LorentzCurve(double x, double t, double s)
    {
        return LorentzCurveFunction(x, t, s);
    }

    /// <summary>
    /// The Neilscheparabola function.
    /// </summary>
    /// <param name="x">The x value for the function.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double NeilscheParabola(double x)
    {
        return (x*Math.Sqrt(x));
    }

    /// <summary>
    /// The Normal distribution function (Gauss).
    /// </summary>
    /// <param name="x">The x value for the function.</param>
    /// <param name="o">The standard degression o.</param>
    /// <param name="m">The expected value m.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double NormalDistributionFunction(double x, double o, double m)
    {
        return Math.Pow((1/(o*Math.Sqrt(2*Math.PI))), (-0.5*Math.Pow(((x - m)/o), 2)));
    }

    /// <summary>
    /// The Normal distribution function (Gauss). Alias command for <see cref="NormalDistributionFunction"/>.
    /// </summary>
    /// <param name="x">The x value for the function.</param>
    /// <param name="o">The standard degression o.</param>
    /// <param name="m">The expected value m.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double NormalDistribution(double x, double o, double m)
    {
        return NormalDistributionFunction(x, o, m);
    }

    /// <summary>
    /// The density function of the standard normal distribution.
    /// </summary>
    /// <param name="x">The x value for the function.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double NormalDistributionDensityFunction(double x)
    {
        return (0.5*Math.PI)*Math.Pow(Math.E, -0.5*(x*x));
    }

    /// <summary>
    /// The density function of the standard normal distribution. Alias command for <see cref="NormalDistributionDensityFunction"/>.
    /// </summary>
    /// <param name="x">The x value for the function.</param>
    /// <returns>The value f(x) at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double NormalDistributionDensity(double x)
    {
        return NormalDistributionDensityFunction(x);
    }

    /// <summary>
    /// The Zeta function.
    /// </summary>
    /// <param name="s">The exponent s for the zetafunction.</param>
    /// <returns>The value for Z(s).</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double ZetaFunction(double s)
    {
        double result = 1;
        double tempuri = 0;
        const double limit = UserDefinedFunctions.StandardIntegrationValue;

        for (int i = 2; i < limit; i++)
        {
            tempuri = 1/(Math.Pow(i, s));
            result += tempuri;

            if (tempuri < limit)
            {
                break;
            }
        }

        return result;
    }

    /// <summary>
    /// The Zeta function. Alias command for <see cref="ZetaFunction"/>.
    /// </summary>
    /// <param name="s">The exponent s for the zetafunction.</param>
    /// <returns>The value for Z(s).</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double Zeta(double s)
    {
        return ZetaFunction(s);
    }

    /// <summary>
    /// The Trigamma function.
    /// </summary>
    /// <param name="z">The z value for the function.</param>
    /// <returns>The value f(z) at the specified position z.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double TrigammaFunction(double z)
    {
        double result = 0;
        const double limit = UserDefinedFunctions.StandardIntegrationValue;

        for (int n = 1; n <= limit; n++)
        {
            result += (1/Math.Pow(z + n, 2));
        }

        return result;
    }

    /// <summary>
    /// The Trigamma function. Alias command for <see cref="TrigammaFunction"/>.
    /// </summary>
    /// <param name="z">The z value for the function.</param>
    /// <returns>The value f(z) at the specified position z.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double Trigamma(double z)
    {
        return TrigammaFunction(z);
    }

    /// <summary>
    /// Computes a approximation result of the gamma function by the Stirling formula for 
    /// values x > 0.
    /// </summary>
    /// <param name="x">The specified number x, for which the gamma function should be calculated.</param>
    /// <returns>Returns the approximation value of the gamma function for the specified number. </returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double GammaFunctionApproximation(double x)
    {
        return Math.Sqrt(2*Math.PI)*Math.Pow(x, x - 0.5)*Math.Pow(Math.E, -x);
    }

    /// <summary>
    /// Computes a approximation result of the gamma function by the Stirling formula for 
    /// values x > 0. Alias command for <see cref="GammaFunctionApproximation"/>.
    /// </summary>
    /// <param name="x">The specified number x, for which the gamma function should be calculated.</param>
    /// <returns>Returns the approximation value of the gamma function for the specified number. </returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double GammaApprox(double x)
    {
        return GammaFunctionApproximation(x);
    }

    /// <summary>
    /// The Gamma function for real numbers.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <returns>The result of the real gamma function at position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double GammaFunction(double x)
    {
        return SpecialFunctionHelper.GammaFunctionHelper(x);
    }

    /// <summary>
    /// The Gamma function for real numbers. Alias command for <see cref="GammaFunction"/>.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <returns>The result of the real gamma function at position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double Gamma(double x)
    {
        return SpecialFunctionHelper.GammaFunctionHelper(x);
    }

    /// <summary>
    /// The lower incomplete Gammafunction T(a,x).
    /// </summary>
    /// <param name="a">The value for evaluation a.</param>
    /// <param name="x">The value for evaluation x.</param>
    /// <returns>
    /// The value ln T(a,x) at the specified position.
    /// </returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double LowerIncompleteGammaFunction(double a, double x)
    {
        double result = 0;
        double tempuri1 = 0;
        double tempuri2 = 0;

        SpecialFunctionHelper.IncompleteGammaFunctionHelper(a, x, ref result, ref tempuri1, ref tempuri2);

        return result;
    }

    /// <summary>
    /// The lower incomplete Gammafunction T(a,x). Alias command for <see cref="LowerIncompleteGammaFunction"/>.
    /// </summary>
    /// <param name="a">The value for evaluation a.</param>
    /// <param name="x">The value for evaluation x.</param>
    /// <returns>
    /// The value ln T(a,x) at the specified position.
    /// </returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double LowerIncompleteGamma(double a, double x)
    {
        return LowerIncompleteGammaFunction(a, x);
    }

    /// <summary>
    /// The upper incomplete Gammafunction T(a,x).
    /// </summary>
    /// <param name="a">The value for evaluation a.</param>
    /// <param name="x">The value for evaluation x.</param>
    /// <returns>
    /// The value ln T(a,x) at the specified position.
    /// </returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double UpperIncompleteGammaFunction(double a, double x)
    {
        double result = 0;
        double tempuri1 = 0;
        double tempuri2 = 0;

        SpecialFunctionHelper.IncompleteGammaFunctionHelper(a, x, ref tempuri1, ref tempuri2, ref result);

        return result;
    }

    /// <summary>
    /// The upper incomplete Gammafunction T(a,x). Alias command for <see cref="UpperIncompleteGammaFunction"/>.
    /// </summary>
    /// <param name="a">The value for evaluation a.</param>
    /// <param name="x">The value for evaluation x.</param>
    /// <returns>
    /// The value ln T(a,x) at the specified position.
    /// </returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double UpperIncompleteGamma(double a, double x)
    {
        return UpperIncompleteGammaFunction(a, x);
    }

    /// <summary>
    /// The incomplete Gammafunction T(a,x).
    /// </summary>
    /// <param name="a">The value for evaluation a.</param>
    /// <param name="x">The value for evaluation x.</param>
    /// <returns>
    /// The value ln T(a,x) at the specified position.
    /// </returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double IncompleteGammaFunction(double a, double x)
    {
        double result = 0;

        if ((a < 0.0) || (x < 0))
        {
            return result;
        }

        double xam = -x + a*Math.Log(x);

        if ((xam > 700) || (a > 170.0))
        {
            return result;
        }

        if (x == 0.0)
        {
            return result;
        }

        if (x <= 1.0 + a)
        {
            double s = 1.0/a;
            double r = s;

            for (int k = 1; k <= 60; k++)
            {
                r *= x/(a + k);
                s += r;

                if (Math.Abs(r/s) < 1e-15)
                {
                    break;
                }
            }

            result = (Math.Exp(xam)*s)/GammaFunction(a);
        }
        else
        {
            double t0 = 0.0;

            for (int k = 60; k >= 1; k--)
            {
                t0 = (k - a)/(1.0 + k/(x + t0));
            }

            result = 1.0 - (Math.Exp(xam)/(x + t0))/GammaFunction(a);
        }

        return result;
    }

    /// <summary>
    /// The incomplete Gammafunction T(a,x). Alias command for <see cref="IncompleteGammaFunction"/>.
    /// </summary>
    /// <param name="a">The value for evaluation a.</param>
    /// <param name="x">The value for evaluation x.</param>
    /// <returns>
    /// The value ln T(a,x) at the specified position.
    /// </returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double IncompleteGamma(double a, double x)
    {
        return IncompleteGammaFunction(a, x);
    }

    /// <summary>
    /// The Loggamma function ln T(x).
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <returns>The value ln T(x) at the specified position.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double LogGammaFunction(double x)
    {
        //Based on the algorithm of C. Bond. Copyright (c) 2003, C. Bond. All rights reserved.
        int n = 0;
        double x0 = x;
        double[] a = new double[]
                         {
                             8.333333333333333e-02,
                             -2.777777777777778e-03,
                             7.936507936507937e-04,
                             -5.952380952380952e-04,
                             8.417508417508418e-04,
                             -1.917526917526918e-03,
                             6.410256410256410e-03,
                             -2.955065359477124e-02,
                             1.796443723688307e-01,
                             -1.39243221690590
                         };

        if (x <= 0.0)
        {
            return 1e308;
        }
        if ((x == 1.0) || (x == 2.0))
        {
            return 0.0;
        }
        if (x <= 7.0)
        {
            n = (int) (7 - x);
            x0 = x + n;
        }

        double x2 = 1.0/(x0*x0);
        double gl0 = a[9];
        const double xp = 2.0*Math.PI;

        for (int k = 8; k >= 0; k--)
        {
            gl0 = gl0*x2 + a[k];
        }

        double gl = gl0/x0 + 0.5*Math.Log(xp) + (x0 - 0.5)*Math.Log(x0) - x0;

        if (x <= 7.0)
        {
            for (int k = 1; k <= n; k++)
            {
                gl -= Math.Log(x0 - 1.0);
                x0 -= 1.0;
            }
        }

        return gl;
    }

    /// <summary>
    /// The Loggamma function ln T(x). Alias command for <see cref="LogGammaFunction"/>.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <returns>The value ln T(x) at the specified position.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double LogGamma(double x)
    {
        return LogGammaFunction(x);
    }

    /// <summary>
    /// The Beta function.
    /// </summary>
    /// <param name="x">The first parameter x of the beta function.</param>
    /// <param name="y">The second parameter y of the beta function.</param>
    /// <returns>The result B(x,y) of the beta function.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double BetaFunction(double x, double y)
    {
        return (GammaFunction(x)*GammaFunction(y))/
               GammaFunction(x + y);
    }

    /// <summary>
    /// The Beta function. Alias command for <see cref="BetaFunction"/>.
    /// </summary>
    /// <param name="x">The first parameter x of the beta function.</param>
    /// <param name="y">The second parameter y of the beta function.</param>
    /// <returns>The result B(x,y) of the beta function.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double Beta(double x, double y)
    {
        return BetaFunction(x, y);
    }

    /// <summary>
    /// The regularized incomplete beta function.
    /// </summary>
    /// <param name="x">The first paramter of the function.</param>
    /// <param name="a">The second paramter of the function.</param>
    /// <param name="b">The tird paramter of the function.</param>
    /// <returns>The value at the specified position x.</returns>
    public static double RegularizedIncompleteBetaFunction(double x, double a, double b)
    {
        if (x == 0.0 || x == 1.0)
        {
            return x;
        }
        else
        {
            bool recur;
            double f;

            if (x > 0.5)
            {
                f = a;
                a = b;
                b = f;
                x = 1.0 - x;
                recur = true;
            }
            else
            {
                recur = false;
            }

            int n = 1;
            int m = 0;
            double g = 0;
            double fn1 = 1;
            double fn2 = 0;
            double gn1 = 1;
            double gn2 = 1;
            double pq = a + b;
            bool neven = false;

            f = 1;

            do
            {
                double dn;

                if (neven)
                {
                    m++;
                    dn = m*x*(b - m)/(a + n - 1.0)/(a + n);
                }
                else
                {
                    dn = -x*(a + m)*(pq + m)/(a + n - 1.0)/(a + n);
                }

                g = f;

                double fn = fn1 + dn*fn2;
                double gn = gn1 + dn*gn2;

                neven = (!neven);
                f = fn/gn;
                fn2 = fn1;
                fn1 = fn;
                gn2 = gn1;
                gn1 = gn;
                n++;
            } while (Math.Abs((f - g)/f) > 1e-10);

            f = f*Math.Pow(x, a)*Math.Pow(1.0 - x, b)*GammaFunction(a + b)/GammaFunction(a + 1.0)/GammaFunction(b);

            if (recur)
            {
                f = 1.0 - f;
            }

            return f;
        }
    }

    /// <summary>
    /// The regularized incomplete beta function. Alias command for <see cref="RegularizedIncompleteBetaFunction"/>.
    /// </summary>
    /// <param name="x">The first paramter of the function.</param>
    /// <param name="a">The second paramter of the function.</param>
    /// <param name="b">The tird paramter of the function.</param>
    /// <returns>The value at the specified position x.</returns>
    public static double RegularizedIncompleteBeta(double x, double a, double b)
    {
        return RegularizedIncompleteBetaFunction(x, a, b);
    }

    /// <summary>
    /// The Central Beta function.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <returns>The result of Central Beta function at position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double CentralBetaFunction(double x)
    {
        return BetaFunction(x, x);
    }

    /// <summary>
    /// The Central Beta function. Alias command for <see cref="CentralBetaFunction"/>.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <returns>The result of Central Beta function at position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double CentralBeta(double x)
    {
        return CentralBetaFunction(x);
    }

    /// <summary>
    /// Computes a approximation result of the Beta function by the Stirling formula for 
    /// values x > 0.
    /// </summary>
    /// <param name="x">The first parameter x of the beta function.</param>
    /// <param name="y">The second parameter y of the beta function.</param>
    /// <returns>The approximated result B(x,y) of the beta function.</returns>
    public static double BetaFunctionApproximation(double x, double y)
    {
        return Math.Sqrt(2*Math.PI)*
               ((Math.Pow(x, x - 0.5)*Math.Pow(y, y - 0.5))/(Math.Pow(x + y, x + y - 0.5)));
    }

    /// <summary>
    /// Computes a approximation result of the Beta function by the Stirling formula for 
    /// values x > 0. Alias command for <see cref="BetaFunctionApproximation"/>.
    /// </summary>
    /// <param name="x">The first parameter x of the beta function.</param>
    /// <param name="y">The second parameter y of the beta function.</param>
    /// <returns>The approximated result B(x,y) of the beta function.</returns>
    public static double BetaFunctionApprox(double x, double y)
    {
        return BetaFunctionApproximation(x, y);
    }

    /// <summary>
    /// The Dirichlet beta function.
    /// </summary>
    /// <param name="s">The value for evaluation.</param>
    /// <returns>The result of the Dirichlet beta function at position s.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double DirichletBetaFunction(double s)
    {
        double result = 0;

        for (int n = 0; n < UserDefinedFunctions.StandardIntegrationValue; n++)
        {
            result += (Math.Pow(-1, n))/(Math.Pow(2.0*n + 1.0, s));
        }

        return result;
    }

    /// <summary>
    /// The Dirichlet beta function. Alias command for <see cref="DirichletBetaFunction"/>.
    /// </summary>
    /// <param name="s">The value for evaluation.</param>
    /// <returns>The result of the Dirichlet beta function at position s.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double DirichletBeta(double s)
    {
        return DirichletBetaFunction(s);
    }

    /// <summary>
    /// The Hurwitz Zeta function, which is define for s > 1.
    /// </summary>
    /// <param name="s">The first parameter s of the function.</param>
    /// <param name="a">The second parameter a of the function.</param>
    /// <param name="limit">The limit of integration.</param>
    /// <returns>The result f(s,a) of the Hurwitz zeta function.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double HurwitzZetaFunction(double s, double a, int limit)
    {
        double result = 0;

        for (int k = 0; k < limit; k++)
        {
            result += 1.0/Math.Pow(k + a, s);
        }

        return result;
    }

    /// <summary>
    /// The Hurwitz Zeta function, which is define for s > 1. Alias command for <see cref="HurwitzZetaFunction"/>.
    /// </summary>
    /// <param name="s">The first parameter s of the function.</param>
    /// <param name="a">The second parameter a of the function.</param>
    /// <param name="limit">The limit of integration.</param>
    /// <returns>The result f(s,a) of the Hurwitz zeta function.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double HurwitzZeta(double s, double a, int limit)
    {
        return HurwitzZetaFunction(s, a, limit);
    }

    /// <summary>
    /// The Lerch Transcendent.
    /// </summary>
    /// <param name="z">The first parameter z.</param>
    /// <param name="s">The second parameter s.</param>
    /// <param name="a">The third parameter a.</param>
    /// <param name="limit">The limit of integration.</param>
    /// <returns>The result f(z,s,a) of the Lerch Transcendent.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double LerchTranscendent(double z, double s, double a, int limit)
    {
        double result = 0;

        for (int k = 0; k < limit; k++)
        {
            result += Math.Pow(z, k)/Math.Pow(k + a, s);
        }

        return result;
    }

    /// <summary>
    /// The Polygamma function.
    /// </summary>
    /// <param name="z">The first parameter z of the polygamma function.</param>
    /// <param name="n">The second parameter n of the polygamma function.</param>
    /// <param name="limit">The limit of integration.</param>
    /// <returns>The result fn(z) of the Polygamma function.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double PolyGammaFunction(double z, int n, int limit)
    {
        double result = 0;

        for (int k = 0; k < limit; k++)
        {
            result += 1.0/Math.Pow(z + k, n + 1);
        }

        return Math.Pow(-1, n + 1)*Factorials.Factorial((uint) n)*result;
    }

    /// <summary>
    /// The Polygamma function. Alias command for <see cref="PolyGammaFunction"/>.
    /// </summary>
    /// <param name="z">The first parameter z of the polygamma function.</param>
    /// <param name="n">The second parameter n of the polygamma function.</param>
    /// <param name="limit">The limit of integration.</param>
    /// <returns>The result fn(z) of the Polygamma function.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double PolyGamma(double z, int n, int limit)
    {
        return PolyGammaFunction(z, n, limit);
    }

    /// <summary>
    /// The divisor function.
    /// </summary>
    /// <param name="n">The specified value n for evaluation.</param>
    /// <param name="x">The xth powers of the divisors of n.</param>
    /// <returns>The result of the divisor function ox(n).</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static int DivisorFunction(int n, int x)
    {
        double result = 0;

        for (uint u = 1; u <= n; u++)
        {
            double tempuri = 0;

            for (int v = 1; v <= u; v++)
            {
                tempuri += Math.Cos((2*Math.PI*v*n)/u);
            }

            result += Math.Pow(u, x - 1)*tempuri;
        }

        return (int) Math.Ceiling(result);
    }

    /// <summary>
    /// The Sigma function.
    /// </summary>
    /// <param name="n">The specified value n for evaluation.</param>
    /// <returns>The result of the Sigma function o1(n).</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static int SigmaFunction(int n)
    {
        return UserDefinedFunctions.DivisorFunction(n, 1);
    }

    /// <summary>
    /// The Sigma function. Alias command for <see cref="SigmaFunction"/>.
    /// </summary>
    /// <param name="n">The specified value n for evaluation.</param>
    /// <returns>The result of the Sigma function o1(n).</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static int Sigma(int n)
    {
        return SigmaFunction(n);
    }

    /// <summary>
    /// The Errorfunction. 
    /// </summary>
    /// <param name="z">The value for evaluation.</param>
    /// <returns>The result of the function at the specified position z.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double ErrorFunction(double z)
    {
        double t = 1.0/(1.0 + 0.47047*Math.Abs(z));
        double poly = t*(0.3480242 + t*(-0.0958798 + t*(0.7478556)));
        double ans = 1.0 - poly*Math.Exp(-z*z);

        if (z >= 0)
        {
            return ans;
        }

        return -ans;
    }

    /// <summary>
    /// The Errorfunction. Alias command for <see cref="ErrorFunction"/>.
    /// </summary>
    /// <param name="z">The value for evaluation.</param>
    /// <returns>The result of the function at the specified position z.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double Error(double z)
    {
        return ErrorFunction(z);
    }

    /// <summary>
    /// The Errorfunction. Alias command for <see cref="ErrorFunction"/>.
    /// </summary>
    /// <param name="z">The value for evaluation.</param>
    /// <returns>The result of the function at the specified position z.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double Erf(double z)
    {
        return ErrorFunction(z);
    }

    /// <summary>
    /// The complementary Errorfunction.
    /// </summary>
    /// <param name="z">The value for evaluation.</param>
    /// <returns>The result of the function at the specified position z.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double ComplementaryErrorFunction(double z)
    {
        return 1 - ErrorFunction(z);
    }

    /// <summary>
    /// The complementary Errorfunction. Alias command for <see cref="ComplementaryErrorFunction"/>.
    /// </summary>
    /// <param name="z">The value for evaluation.</param>
    /// <returns>The result of the function at the specified position z.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double ComplementaryError(double z)
    {
        return ComplementaryErrorFunction(z);
    }

    /// <summary>
    /// The complementary Errorfunction. Alias command for <see cref="ComplementaryErrorFunction"/>.
    /// </summary>
    /// <param name="z">The value for evaluation.</param>
    /// <returns>The result of the function at the specified position z.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double Cerf(double z)
    {
        return ComplementaryErrorFunction(z);
    }

    /// <summary>
    /// The WhittakerM function Mk,m(z). Its the Whittaker function of the first kind.
    /// </summary>
    /// <param name="k">The value of k.</param>
    /// <param name="m">The value of m.</param>
    /// <param name="z">The value of z.</param>
    /// <returns>The result of the function at the specified position.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double WhittakerMFunction(int k, int m, double z)
    {
        double result = 0;

        for (uint n = 0; n <= 10; n++)
        {
            result += (PochhammerSymbol(m - k + 0.5, n)
                       /(Factorials.Factorial(n)*PochhammerSymbol(2*m + 1, n)))
                      *Math.Pow(z, n);
        }

        return Math.Pow(z, 0.5 + m)*Math.Pow(Math.E, -z/2)*result;
    }

    /// <summary>
    /// The WhittakerM function Mk,m(z). Its the Whittaker function of the first kind. 
    /// Alias command for <see cref="WhittakerMFunction"/>.
    /// </summary>
    /// <param name="k">The value of k.</param>
    /// <param name="m">The value of m.</param>
    /// <param name="z">The value of z.</param>
    /// <returns>The result of the function at the specified position.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double WhittakerM(int k, int m, double z)
    {
        return WhittakerMFunction(k, m, z);
    }

    /// <summary>
    /// The confluent hypergeometric function F1(a; b; z) of the first kind.
    /// </summary>
    /// <param name="a">The value of a.</param>
    /// <param name="b">The value of b.</param>
    /// <param name="z">The value of z.</param>
    /// <returns>The result of the function at the specified position.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double ConfluentHypergeometricFunctionFirstKind(double a, double b, double z)
    {
        double result = 0;
        UInt64[] factorials = Factorials.FactorialList20();

        for (UInt64 k = 0; k < 20; k++)
        {
            result += (PochhammerSymbol(a, k)/PochhammerSymbol(b, k))
                      *(Math.Pow(z, k)/factorials[k]);
        }

        return result;
    }

    /// <summary>
    /// The confluent hypergeometric function F1(a; b; z) of the first kind. 
    /// Alias command for <see cref="ConfluentHypergeometricFunctionFirstKind"/>.
    /// </summary>
    /// <param name="a">The value of a.</param>
    /// <param name="b">The value of b.</param>
    /// <param name="z">The value of z.</param>
    /// <returns>The result of the function at the specified position.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double ConfluentHypergeometricFirstKind(double a, double b, double z)
    {
        return ConfluentHypergeometricFunctionFirstKind(a, b, z);
    }

    /// <summary>
    /// The Pochhammer symbol (x)n.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <param name="n">The order n for the Pochhammer symbol.</param>
    /// <returns>The result of the function at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double PochhammerSymbol(double x, double n)
    {
        return GammaFunction(x + n)/GammaFunction(x);
    }

    /// <summary>
    /// The Kummerfunction 2F1(1/2 + m - k, -n; 2m + 1; 1) of the first kind.
    /// </summary>
    /// <param name="k">The value of k.</param>
    /// <param name="m">The value of m.</param>
    /// <param name="n">The value of n.</param>
    /// <returns>The result of the function at the specified position.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double KummerFunctionFirstKind(double k, double m, double n)
    {
        double a = GammaFunction(2*m + 1);
        double b = GammaFunction(m + 0.5 + k + n);
        double c = GammaFunction(m + 0.5 + k);
        double d = GammaFunction(2*m + 1 + n);

        return (a*b)/(c*d);
    }

    /// <summary>
    /// The Kummerfunction 2F1(1/2 + m - k, -n; 2m + 1; 1) of the first kind. Alias command for <see cref="KummerFunctionFirstKind"/>.
    /// </summary>
    /// <param name="k">The value of k.</param>
    /// <param name="m">The value of m.</param>
    /// <param name="n">The value of n.</param>
    /// <returns>The result of the function at the specified position.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double KummerFirstKind(double k, double m, double n)
    {
        return KummerFunctionFirstKind(k, m, n);
    }

    /// <summary>
    /// Bessel function Jn(x) of the first kind.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <param name="n">The order n of the function.</param>
    /// <returns>The result of the function at the specified position x.</returns>
    public static double BesselFunctionFirstKind(double x, int n)
    {
        bool isNull = false;

        if (n == 0)
        {
            n++;
            isNull = true;
        }

        double[] jn = new double[n + 1];
        double[] yn = new double[n + 1];
        double[] jnp = new double[n + 1];
        double[] ynp = new double[n + 1];
       
        SpecialFunctionHelper.BesselFunction(n, x, jn, yn, jnp, ynp);

        if (isNull)
        {
            return jn[0];
        }

        return jn[n];
    }

    /// <summary>
    /// Bessel function Jn(x) of the first kind. Alias command for <see cref="BesselFunctionFirstKind"/>.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <param name="n">The order n of the function.</param>
    /// <returns>The result of the function at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double BesselFirstKind(double x, int n)
    {
        return BesselFunctionFirstKind(x, n);
    }

    /// <summary>
    /// Bessel function Yn(x) of the second kind.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <param name="n">The order n of the function.</param>
    /// <returns>The result of the function at the specified position x.</returns>
    public static double BesselFunctionSecondKind(double x, int n)
    {
        bool isNull = false;

        if (n == 0)
        {
            n++;
            isNull = true;
        }

        double[] jn = new double[n + 1];
        double[] yn = new double[n + 1];
        double[] jnp = new double[n + 1];
        double[] ynp = new double[n + 1];

        SpecialFunctionHelper.BesselFunction(n, x, jn, yn, jnp, ynp);

        if (isNull)
        {
            return yn[0];
        }

        return yn[n];
    }

    /// <summary>
    /// Bessel function Yn(x) of the second kind. Alias command for <see cref="BesselFunctionSecondKind"/>.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <param name="n">The order n of the function.</param>
    /// <returns>The result of the function at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double BesselSecondKind(double x, int n)
    {
        return BesselFunctionSecondKind(x, n);
    }

    /// <summary>
    /// The PSI function.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <returns>The result of the function at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double PSIFunction(double x)
    {
        int n;
        double ps, x2;

        double[] a = {
                         -0.8333333333333e-01,
                         0.83333333333333333e-02,
                         -0.39682539682539683e-02,
                         0.41666666666666667e-02,
                         -0.75757575757575758e-02,
                         0.21092796092796093e-01,
                         -0.83333333333333333e-01,
                         0.4432598039215686
                     };

        double xa = Math.Abs(x);
        double s = 0.0;
        const double el = 0.5772156649015329;

        if ((x == (int) x) && (x <= 0.0))
        {
            ps = 1e308;

            return ps;
        }

        if (xa == (int) xa)
        {
            n = (int) xa;

            for (int k = 1; k < n; k++)
            {
                s += 1.0/k;
            }

            ps = s - el;
        }
        else if ((xa + 0.5) == ((int) (xa + 0.5)))
        {
            n = (int) (xa - 0.5);

            for (int k = 1; k <= n; k++)
            {
                s += 1.0/(2.0*k - 1.0);
            }

            ps = 2.0*s - el - 1.386294361119891;
        }
        else
        {
            if (xa < 10.0)
            {
                n = 10 - (int) xa;

                for (int k = 0; k < n; k++)
                {
                    s += 1.0/(xa + k);
                }

                xa += n;
            }

            x2 = 1.0/(xa*xa);
            ps = Math.Log(xa) - 0.5/xa +
                 x2*
                 (((((((a[7]*x2 + a[6])*x2 + a[5])*x2 + a[4])*x2 + a[3])*x2 + a[2])*x2 + a[1])*x2 +
                  a[0]);
            ps -= s;
        }

        if (x < 0.0)
        {
            ps -= Math.PI*Math.Cos(Math.PI*x)/Math.Sin(Math.PI*x) - 1.0/x;
        }

        return ps;
    }

    /// <summary>
    /// The PSI function. Alias command for <see cref="PSIFunction"/>.
    /// </summary>
    /// <param name="x">The value for evaluation.</param>
    /// <returns>The result of the function at the specified position x.</returns>
    [Microsoft.SqlServer.Server.SqlFunction]
    public static double PSI(double x)
    {
        return PSIFunction(x);
    }
} ;